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Abstract 


Phase field model (PFM) is an efficient fracture modeling method and has high potential 
for hydraulic fracturing (HF). However, the current PFMs in HF do not consider well 
the effect of in-situ stress field and the numerical examples of porous media with stress 
boundary conditions were rarely presented. The main reason is that if the remote stress 
is applied on the boundaries of the calculation domain, there will be relatively large de- 
formation induced on these stress boundaries, which is not consistent with the engineering 
observations. To eliminate this limitation, this paper proposes a new phase field method 
to describe quasi-static hydraulic fracture propagation in porous media subjected to stress 
boundary conditions, and the new method is more in line with engineering practice. A 


new energy functional, which considers the effect of initial in-situ stress field, is established 


and then it is used to achieve the governing equations for the displacement and phase 
fields through the variational approach. Biot poroelasticity theory is used to couple the 
fluid pressure field and the displacement field while the phase field is used for determin- 
ing the fluid properties from the intact domain to the fully broken domain. In addition, 
we present several 2D and 3D examples to show the effects of in-situ stress on hydraulic 
fracture propagation. The numerical examples indicate that under stress boundary condi- 
tion our approach obtains correct displacement distribution and it is capable of capturing 


complex hydraulic fracture growth patterns. 
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1 Introduction 


Hydraulic fracture propagation in porous media is one of the most attractive and significant 
research topics in mechanical, geological, energy, and environmental engineering. In particular, 
hydraulic fracturing (HF) [1] has been widely used to exploit oil, tight gas, and shale gas from 
reservoirs that were unexploitable in past decades. The main reason for this is that the injection 
of highly pressurized fluid into a reservoir forms a fracture network for resource transportation. 
Another application of HF is the measurement of in-situ stress [2]. In addition, HF can be 
also applied in an enhanced geothermal system to accelerate heat extraction [3]. Nevertheless, 
despite its advantages, HF still brings some controversies because the fracturing fluid may leak 
and further contaminate the underground space and surface due to unfavorable fracture growth 


paths [4, 5]. Therefore, an accurate numerical tool is critically important in the prediction of 


complex hydraulic fracture propagation in porous media. 

However, correctly modeling hydraulic fracture in porous media is challenging and full of 
complexity due to solid-fluid interaction, fracture network, and different boundary conditions. 
This has prompted the development of various numerical methods for modeling fracturing pro- 
cesses, and these approaches can be classified into two types in the continuum framework: the 
discrete method and the smeared method. The discrete methods introduce displacement. dis- 
continuity for fractures and among the most popular are the extended finite element method 
(XFEM) [6, 7], generalized finite element method (GFEM) [8], boundary element method [9-11], 
phantom-node method [12, 13], and element-erosion method [14, 15]. On the other hand, in the 
smeared methods the displacement is continuous across a fracture and the gradient damage 
model [16], the screened Poisson method [17], and the phase field model (PFM) [18-22] are the 
best-known. 

Among all the fracture modeling methods, the PFMs are attracting more and more attention 
in recent years. In this method, an additional scalar field ¢ € [0, 1] is used to reflect the extent 
of fracture where ¢ = 0 represents an intact material and ¢ = 1 indicates a fully broken material 
(some literatures [23] use a phase field s = 1 — ¢ and s = 0, 1 represents the fully broken state 
and intact state, respectively). In addition, the transition zone with @ € (0,1) is controlled 
by an intrinsic length scale parameter lọ. After being first proposed by Bourdin et al. [18], 
the PFM was further promoted by Bourdin et al. [18], Miehe et al. [19, 20], Borden et al. [24], 
Hofacker and Miehe [25, 26]. Compared with XIGA [27], XFEM [28, 29], cohesive zone model]|30, 
31], and continuum damage model [32], PFM has ease of implementation. Furthermore, the 
recent development has shown that PFMs have a high efficiency in predicting complex fracture 


propagation patterns, even in 3D situations due to these reasons: i) simulations can be performed 


on a fixed mesh without any remeshing or adaptive technique; ii) complex fracture patterns such 
as branching and junction are automatically captured; iii) external fracture criteria or fracture 
surface tracking algorithms are not required; iv) PFMs can easily simulate fracture propagation 
in heterogeneous media, and v) no penetration criteria are required for hydraulic fracture when 
a layer interface is encountered. 

The aforementioned advantages have also contributed to the development of PFMs in hy- 
draulic fracturing. For example, in recent years many researchers have tried to couple the PFMs 
to HF and made some progress [33-45]. However, the presented examples in these contribu- 
tions all only established homogeneous Dirichlet boundary conditions for the displacement field. 
The current PFMs in HF therefore cannot consider well the effect of in-situ stress on fracture 
propagation. The main reason is that if the remote stress is applied on the boundaries of the 
calculation domain, there will be relatively large deformation induced on these stress boundaries, 
which is not consistent with the engineering observations. In fact, in geological environment, 
the displacement on the boundaries where the remote stress acts should be zero before HF. It 
should be noted that although a recently developed PFM [46] attempted to analyze the effect of 
stress boundary condition on fracture pattern, the main drawback remained unsolved because a 
large deformation was observed on the stress boundaries. 

To eliminate the limitation from the initial stress field, this paper proposes a new phase field 
method to describe quasi-static hydraulic fracture propagation in porous media subjected to 
stress boundary conditions. A new energy functional, which fully considers the effect of initial 
in-situ stress field, is established and then it is used to achieve the governing equations of strong 
form for the displacement and phase fields through the variational approach. Biot poroelasticity 


theory is used to couple the fluid pressure field and the displacement field while the phase field 


is used for determining the fluid properties from the intact domain to the fully broken domain. 
In addition, several 2D and 3D examples are presented to show the effects of in-situ stress on 
hydraulic fracture propagation. The numerical examples indicate that under stress boundary 
conditions the new PFM obtains correct displacement distribution and it is capable of capturing 
well complex hydraulic fracture growth patterns. 

This paper is organized as follows. Section 2 describes the mathematical model for the new 
PFM and Section 3 shows the global numerical algorithm using a staggered scheme. In Sections 
4 and 5 the 2D and 3D examples are presented to demonstrate the capability of the new model. 


The present work and outlook for future development are concluded in Section 6. 


2 Mathematical model 


2.1 New energy functional 


Let us consider a poro-elastic domain Q C R? (d € {2,3}) in Fig. 1 where the external and 


internal discontinuity boundaries are denoted as ðQ and I, respectively. The porous medium is 


assumed homogeneous and isotropic with compressible and viscous fluid in the pores. Let T > 0 


be the computational time interval and u(x,t) C R? be the displacement field at time t € [0, T] 


and the position æ. In addition, the displacement field must satisfy the time-dependent Dirichlet 
boundary conditions, u;(a,t) = gi(a,t), on OQ; € OQ, and the Neumann condition on Op. 
The stress, displacement, and fluid pressure of the porous domain are shown in Fig. 2. 
Because of the long-term consolidation or other geological effects, the domain has formed an 
initial stress field oo, initial displacement field uo, and fluid pressure po, as shown in Fig. 2a. 


On the other hand, in Fig. 2b, after the domain is excavated or fractured by fluid injection, 


(a) Sharp fracture (b) Diffusive fracture 


Fig. 1 Sharp and diffusive fractures in the poro-elastic medium 


the stress, displacement, and fluid pressure are o, Up + u, and po + p, respectively. Here, 
the displacement u and pressure p are the induced relative displacement and fluid pressure by 
engineering activities such as HF. In this work, our presented method is only for calculating o, 
u, p, and the hydraulic fracture pattern in the case that og, Uo, and po are known in advance, 
while how to obtain these initial fields is beyond the scope of this paper. 

It should be noted that for the porous medium in the geological environment, the displace- 
ment ug resulting from long-term geo-stress is always ignored in the stability analysis for un- 
derground engineering [47, 48] and only the displacement u caused by fracture formation or 
engineering excavation is calculated. That is, the porous medium Q is subjected to an initial 
stress field op» and if the stress in Q is equal to oo, the displacement field must be 0. In addition, 
the initial fluid pressure po is set as 0 for the purpose of simplicity and because the fluid pressure 
p results only from the relative displacement u and fluid injection [45]. 

If the body force is ignored, the basic idea of the previous phase field models for porous- 


elastic media [45] is to construct an energy functional Ų composed of the elastic energy V.(e), 


Or 


(a) Initial state (b) After excavation or HF 


Fig. 2 Stress, displacement and fluid pressure in geological environment 


fracture energy Wy, external work West, and the energy contribution of fluid pressure p: 


V(u,p,T) = w.(e)dQ — I ap- (V -u)dQ + f Gdr — fi- uds (1) 


AL Ot 


VA pressure-related term Wy West 


where w.-(€) is the elastic energy density; a € (0,1] is the Biot coefficient; Ge is the critical 
energy release rate; f; is the traction on the Neumann boundary on dS. In addition, the linear 


strain tensor € is given by 


and for an isotropic linear elastic medium, the elastic energy density %e(€) reads [20] 


1 
pele) = 5 EINES + HEijEij (3) 


where A, u > 0 are the Lamé constants and 


Ev 
N= (1+ v)(1 — 2v) (4) 
= E 
K= Frag) 


with E and v being Young’s modulus and Poisson’s ratio, respectively. 

However, as seen in Eq.(1), the previously used energy functional cannot account for the 
effect of the initial stress field and it is therefore not suitable for modeling fracture propagation 
in a geomaterial under stress boundary condition. An evidence for this can be observed in 
Shiozawa et al. [46] where a stress boundary has a relatively large deformation when the fluid 
pressure is 0. Hence, to deal with the stress boundary and account for the effect of the initial 


stress field, we establish a new energy functional L for the phase field modeling: 


L(u,p, T) = We(e)dQ + | Oo : edQ — I ap: (V- wdo+ f Gar — fi- udS (5) 
Ar Q T an 


r 


Ye Initial stress induced pressure-related term Vy West 


where the sum of the first and second terms means the incremental elastic energy from the state 
of uo to the state of up + u. It will be shown in the following sections that this incremental 
elastic energy drives the fracture propagation in the porous medium Q. In addition, it should 
be noted that by using the new energy functional (5), the fracture deflection phenomenon due 
to stress contrast can be also well captured by the phase field model, as shown in the numerical 


examples presented in this paper. 


or 


Ld 


2.2 Phase field description 


It can be seen from Eqs. (1) and (5) that the energy functional contains a sharp internal surface 
I’, which increases the difficulty in minimizing the energy functional when the variational method 
is used. Therefore, to simplify the numerical implementation, a phase field ¢(a, t) € [0,1] is used 
to smear the sharp fracture as shown in Fig. 1b where ¢ = 0 and ¢ = 1 represent an intact 
material and a fully broken material, respectively. 


For a fracture in a 1D bar, the solution for the phase field is an inverse exponential function 


(a) [20]: 


s) = exp (- E) (6) 


where x = a is the fracture location and lọ denotes the intrinsic length scale parameter. The 
length scale lo is also required for 2D and 3D problems where the crack surface density per unit 


volume is used in terms of the phase field and its gradient [20]: 


$ 


ETE £42 2yo. vg (7) 


Note that the length scale parameter controls the width of a diffused fracture and a larger 
lo shows a lower nominal tensile strength in the phase field modeling [49]. Note that the sharp 
fracture can be recovered if lọ tends to zero, which is known as T-convergence [20]. In addition, 
the length scale is assumed much larger than the pore size of the domain Q in this study and 


therefore the fracture energy in Eq. (5) is rewritten as 


v= [Gare f coan= f a, (= + Siver) an (8) 
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The elastic energy Y, must be decomposed into tensile and compressive parts to ensure 
cracks only under tension [24]. Therefore, we follow the strain spectral decomposition of Miehe 


et al. [20]: 


E4 => (Ea)iMa Q Na (9) 


where e4 are the tensile and compressive strain tensors, respectively; €a and na are the principal 


. are defined as [20]: (*)4 = (* | x |)/2. 


strain and its direction; the operators (*} 


Applying the decomposed strain tensor, the tensile and compressive parts of the elastic energy 


density are written as 


yè le) = S(tr(e))3 + utr (e) (10) 


We follow Borden et al. [24] and the compressive part of the elastic energy density is assumed 


not to affect the fracture propagation. Therefore, the elastic energy is rewritten as 


Y, = Veled = | [9O e) + Yr (e)] 40 (11) 


Qr 


where g(@) is a degradation function, and g(0) = 1, g(1) = 0, g'(1) = 0 must be satisfied [23]. 
Although there are many forms for g(@), a quadratic form of g(¢) = (1—k)(1—¢)? + k is applied 


in this work with k = 107° being a stability parameter to prevent numerical singularity when 
o=0. 


The degradation function is also applied to the energy variation due to the initial stress field: 


f oo: EdQ & f g(o)oo : edQ (12) 
Ar Q 
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which means the initial stress field does not contribute to the energy functional in a fully broken 


region with ¢ = 0. 


2.3 Governing equations for the displacement 


Now, substituting the fracture energy (8), the elastic energy (11), and Eq. (12) into Eq. (5), 


the energy functional can be rewritten as 


p= [ [oot (e) + y7 (e)] a2 + [ g(d)oo : ed — f ap- (V - u)dQ+ 


f G p 2 “V6 ve] dQ— | f,-udS (13) 
Q 


2lo IQ}, 


We then apply the variational approach [50] where fracture initiation and propagation at 
time t € [0,7] is a process to minimize the energy functional L. Therefore, the first variation of 


the energy functional L is set as 0, yielding 


ôL = i (os F g(0) Ooi; = apõi;) Mj — fu] uid S -f (of, + g(0) Ooi; — apõij) du;dQ — 
Qh, 
s.a —) Oren” 


© @ 
j Ge o? ð 
f E (OYT + ooij€ig) + E = Clog saas f (Em) ôpdS =0 (14) 
ee en, 
© ® 


where g/(¢) = dg(¢)/d@ = 2(¢ — 1)(1 — k) and m; is the component of the outward-pointing 


normal vector of the boundary. of; is the component of the effective stress tensor o(€)° induced 
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by the displacement u: 


OVE | Abe 
-} 


= |(1—k)(1—¢)? + k] [A(tr(e)) 4 + 2ue+] + A(tr(e))_I + 2pe- 


o° = g(¢) 


Rexd 


where I is the identity tensor € | 


Now, we define the total stress tensor o as 


a(€)=a%(e)+g9(d)oo-—apI, inQ x (0,7) (16) 


Combining Eqs. (15) and (16), it can be seen that in a pure tension state without fluid, if 
the phase field ¢ = 1, the total stress ø is 0 in the fully broken region. In addition, because 


arbitrary admissible ĝu must satisfy Eq. (14), Eq. (14)@ gives rise to the governing equation: 


00%; 


n 0, inQ x (0,T] (17) 


subjected to the Dirichlet boundary and the Neumann boundary given by 


OijMj = Ses on OO), x (0, T] (18) 


which can be derived from Eq. (14). 
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2.4 Governing equations for the phase field 


Because of the arbitrariness of 6¢, Eq. (14)@ results in the original governing equation for the 


phase field: 


2lo(1 — k) (Y + ToijEiz) bal ae po? = 2lo(1 — k) (WF + doijEig) in Q x (0, T] (19) 


Ge ° Oar? Ge 


However, Eq. (19) cannot ensure the irreversibility condition I(a,s) € Tr(æ,t)(s < t). 
Therefore, a history reference field H is constructed to form a monotonically increasing phase 
field: 


H(a,t) = max |Y} (e(x,s)) + oo: e(a,s)]}, inQ x (0,T] (20) 


s€[0,t] 


Replacing WF + o0:;€i; by H(a,t) in Eq. (19), the strong form for the phase field is given by 


2lo(1 —k)H Od 2(1-k)H. 
ac + 1 o-h Ou? = A , inQ~x (0,T] (21) 
which is subjected to the Neumann condition: 
OF i =0, on ðQ x (0,T] (22) 


Therefore, the evolution equation of the phase field (21) indicates that the fracture initiation 
and propagation is driven by the history reference field H (x,t), which represents the historic 


high incremental elastic energy from the uo state to the uo + u state during the period (0, t]. 
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2.5 Fluid flow in the porous medium 


As described in Subsection 2.1, the fluid in poro-elastic medium is assumed compressible and 


viscous. Then, we calculate the fluid flow in three subdomains, namely the unbroken domain 


(reservoir domain) (,(t), fractured domain Q(t) and transition domain Q(t), respectively. 


These subdomains are distinguished according to the phase field in Table 1 where cı and cy are 


two phase field thresholds. 


Table 1: Subdomain definition 


Subdomain 

Reservoir domain 
Transition domain 
Fractured domain 


Phase field 
esa 
Cı < O < C2 
Q> C 


In this study, we calculate the hydraulic parameters in the transition domain by using a linear 


interpolation from the reservoir and fractured domains. Therefore, two indicator functions X%r 


and x, [46] are established. y, = 1 in the reservoir domain and x, = 0 in the fractured domain. 


In contrast, xf = 0 in the reservoir domain and yf = 1 in the fractured domain. In addition, 


the indicator functions satisfy the following: 


Xr(-, Q) = 2 ? 
2— C4 


Ci << 


CO << 


Darcy’s law is used to model the fluid flow and the mass conservation equation in Zhou et al. 


[45] is applied to the whole calculation domain: 


o 
pS +V + (pv) = dm — pax, 


ot 


14 


Evol 


Ot (24) 


where p, S, V, Evol = V-u, and qm represent the fluid density, storage coefficient, flow velocity, 
volumetric strain, and fluid source term, respectively. In Eq. (24), p = prxr + pfXf and 
a= QrXr + afXf Where p, and py are the fluid densities in the reservoir and fracture domains 
while a, and ay are the Biot coefficients of the reservoir and fractured domains. Naturally, 
af = 1 is set for the fractured domain and therefore a = a,x, + Xf- 

In Eq. (24), the storage coefficient S is expressed as [45] 


(a = Ep)(1 — a) 


S = EpC + Ke 


(25) 


where £p, c, and Ky, are the porosity, fluid compressibility, and bulk modulus of the calcula- 
tion domain, respectively. Denoting c, and cy as the fluid compressibility in the reservoir and 
fractured domains, we have c = c,X, + CFXF- 


The gravity is not considered in this study; therefore the Darcy’s velocity v is calculated as 
v =- Vp (26) 


where Hefs denotes fluid viscosity. Similarly, erp = HrXr + ufXs With ur and us being the fluid 
viscosity in the reservoir and fractured domains, respectively. Keşp is the effective permeability 
and Keff = krXr +kpxp where k, and kp are the permeabilities of the reservoir and fractured 
domains. Replacing Eq. (26) into Eq. (24), the governing equation for fluid flow is expressed in 
terms of the fluid pressure p: 

g _ wy. PRett Evo (27) 


Vp = m r 
uas O= Om = pO 
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which is subjected to the Dirichlet condition on Np and Neumann condition on OQ yn [45]: 


p = pp on üp 


—n-pv= My on Ny 


where pp and My are the prescribed pressure and mass flux. 


3 Numerical algorithm 


We implement the proposed phase field model in the framework of finite element method (FEM) 
and the FE discretization can be seen in Appendix. In addition, a finite difference method is 
applied for the time discretization. The overall procedure of our numerical algorithm for solving 
the three fields is shown in Fig. 3 where a staggered scheme is used. That is, the three fields 


are solved sequentially and independently in each time step. 


Solve the discrete equations of displacement and pressure for (u, p" 
Renew the history reference field for Hg: 


Solve the discrete equation of phase field for g’*! 


Next time step 


Do all the three fields satisfy the tolerance requirement? 


Fig. 3 Numerical scheme in each time step 


To reduce the computational effort, we implement the numerical method in the commercial 
software, COMSOL Multiphysics, which has been proven to have high capability for multi-field 
problems. In all simulations, the maximum iteration number is set as 150 due to high nonlinearity 
resulting from the derivative 0o°/0e and the phase field which degrades the stiffness matrix of 
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the displacement field. In addition, a stabilization and convergence acceleration method—the 
Anderson acceleration technique is applied. For more detail on the COMSOL implementation, 


the readers can be referred to the previous study [51]. 


4 2D examples 


In this section, 2D examples of specimens subjected to internal fluid injection are presented 
to prove the capability of the proposed phase field method. All the pre-existing cracks are 
established by introducing an initial history field with a relatively high H [24] and the source 


term in the pre-existing notches is set as qr = 10 kg/(m?-s) in all examples. 


4.1 Fractures from a horizontal notch 


The first example tests the hydraulic fracture propagation from a horizontal notch in a square 
domain. The geometry and boundary conditions of this example are shown in Fig. 4 with 
increasing fluid volume in the initial notch. All the outer boundaries of the calculation domain 
has a fluid pressure of p = 0 and a zero tangential displacement in order to remove the effect of 
rigid body displacement. In addition, the parameters for calculation are listed in Table 2. 

We discretize the calculation domain with unstructured triangular elements with the maxi- 
mum element size h = 0.05 m. Therefore, the linear shape functions are used for all the three 
physical fields. In addition, the time step At = 0.05s is adopted for the simulation. The initial 
stress state of oro = 0.5 MPa, oy) = 0.5 MPa and 0.9 = v (ozo + Oyo) is applied with v being 
Poisson’s ratio, which can be evaluated through the well-known elasticity relationship from A and 


u. As a comparison, we also model the fracture propagation by using the previously developed 
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Fig. 4 Geometry and boundary conditions of the calculation domain 


Table 2: Basic calculation parameters 


Parameter Value Unit Parameter Value Unit 

u 23.08 GPa À 34.62 GPa 

Ge 500 N/m k 1x 107° - 

lo 0.1 m Cl 0.4 - 

C9 1.0 = EpR 0.05 — 

PR; PF 1.0x 10% kg/m? QR 0.05 = 

qR 0 kg/(m*-s) || qF 0 kg/(m?-s) 
kr 1x107 m? kp 8.333 x 1074 m? 

CR 1x107 1/Pa CF 1 x 1078 1/Pa 

UR 1x107 Pas UF 1 x 1078 Pas 
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method [45] without considering the effect of initial stress field. 

Figures 5 and 6 show the fracture propagation obtained by using the proposed phase field 
method and that without considering the effect of initial stress field. The figures indicate that 
discarding the effect of initial stress does not affect the fracture pattern for the first example; 
in both situation the fracture propagates along the horizontal direction. However, the time for 
fracture initiation and propagation is different and the fracture length is smaller if the initial 


stress field is considered in the constitutive model. 


Fracture half-length = 0.4162 m Fracture half-length = 0.8918 m Fracture half-length = 1.789 m 


(a) t=2.4s (b) t=3.6s (c) t= 4.88s 


Fig. 5 Fracture propagation pattern by using the proposed PFM 


Fracture half-length = 0.4358 m Fracture half-length = 1.389 m Fracture half-length = 2.023 m 


(a) t=2.4s (b) t=3.6s (c) t= 4.88s 


Fig. 6 Fracture propagation pattern without considering the effect of initial stress field 
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The displacement field at time t = 0 s is shown in Fig. 7 where the region of @ > 0.95 is 
removed to reflect the shape of the fully broken domain. It can be observed that the proposed 
PFM and the method of Zhou et al. [45] achieve different initial displacement field for the 
problem of a poro-elastic domain subjected to stress boundary condition. Fig. 7a indicates that 
only a small displacement appear around the center of the notch for the proposed PFM because 
the “excavation” or stiffness degradation of the initial notch produces displacement towards the 
broken domain. However, if the PFM does not consider the effect of initial stress field, all the 
outer boundaries have rather large displacements as shown in Fig. Fig. 7b. In addition, Fig. 8 
compares the displacements along the top and left boundaries at time t = 0 s. As observed, for 
the method of Zhou et al. [45], the stress boundary condition produces large initial displacements 
along the left and top boundaries while the initial displacement on the boundaries is negligible by 
using the proposed PFM. In summary, comparisons in Figs. 7 and 8 indicate the proposed PFM 
will achieve better displacement distribution compared with those PFMs without considering the 
effect of initial stress field, especially for the poro-elastic medium subjected to stress boundary 


condition in the geological environment. 


4 482x10 4 209x107 
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Fig. 7 Comparison of the displacement field at t = 0 s for the proposed PFM and Zhou et al. 
[45] (Unit: m) 
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Fig. 8 Comparison of the displacement along the top and left boundaries at t = 0 s for the 
proposed PFM and Zhou et al. [45] 


Figure 9 shows the effect of initial stress field on the fluid pressure-time curve. Note that 
the data at the center of the initial notch are selected. Similar to the fracture pattern, Fig. 
9 indicates the fluid pressure is only slightly affected by the initial stress field. In the current 
example, the time for fracture initiation is reduced if the initial stress field is not considered; 
therefore the fluid pressure-time curve has a earlier drop stage and a lower maximum pressure 


compared with those obtained by the proposed PFM as shown in Fig. 9. 


The proposed PFM 
Zhou et al. (2018) 


Mid-point fluid pressure [MPa] 


t [s] 


Fig. 9 Comparison of the fluid pressure-time curve for the proposed PFM and Zhou et al. [45] 
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The evolution of the displacement field obtained by the proposed PFM is shown in Fig. 10. 
As expected, the maximum displacement occurs in the center of the fracture and it increases 
with the increasing time. This phenomenon is also consistent with those observations in a 
porous medium with fixed displacement boundaries [36, 45, 52], which indirectly reflects the 


applicability of the proposed PFM in this study. 
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Fig. 10 Evolution of the displacement field by using the proposed PFM (Unit: m) 


It is well-known that the hydraulic fracture pattern is highly affected by the stress contrast 
acting on the outer boundaries of the calculation domain. Therefore, in this example, we change 
the ratio of o40/Czo to [0.5, 1,2, 10] with czo = 0.5 MPa unchanged to demonstrate the effect of 
stress contrast. By using the proposed PFM, the calculated fracture paths at time t = 10 s are 
shown in Fig. 11. It can be observed that for o,9/or0 = 0.5, 1, and 2, the fracture from the 
initial notch propagates horizontally and the fracture length decreases as the ratio of a,0/o20 
increases. However, when 09/070 = 10, the fracture deflects and propagates along the direction 
of the maximum in-situ stress Smar, Which is consistent with the engineering observations in 
hydraulic fracturing. In addition, the effect of the ratio of o,9/7,9 on the fluid pressure-time 
curve is shown in Fig. 12. The maximum fluid pressure at the mid-point of the initial notch is 
observed to increase with the increasing oyo/ox0- 
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a) Fy0/Tx0 = 0.5 b) ayo/ex0 = 1 
c) Fyo/ex0 = 2 d) gyo/ex0 = 10 


Fig. 11 Fracture propagation from a horizontal initial notch at time t = 10 s under different 
Oyo / Oxo 


Mid-point fluid pressure [MPa] 


t [s] 


Fig. 12 Effect of the ratio of oyo/oz9 on fluid pressure-time curve for a horizontal initial notch 
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Figure 13 describes the relative incremental mid-point fluid pressure and effective stress at 
fracture tip under different o,o/o,9 in the case of a horizontal initial notch. Note that the 
effective stress along the y direction is accordingly used to show the effect of ayo, and the 
case of dyo/029 = 0.5 is used as the reference. As shown in Fig. 13, when the vertical stress 
0,0 increases linearly, the mid-point fluid pressure and effective stress at the fracture tip also 
increase linearly. This phenomena indirectly verifies the feasibility and practicability of the 
proposed PFM. Furthermore, the incremental fluid pressure is slightly larger than the vertical 


stress variation while the incremental effective stress at fracture tip is slightly lower. 


Mid-point fluid pressure 


— p 


y0 


=- =- Effective stress o at fracture tip 


Incremental pressure or stress [MPa] 


Oa / Oa 


Fig. 13 Relative incremental fluid pressure and effective stress at fracture tip under different 
00/0 in the case of a horizontal initial notch 


4.2 Fracture from an inclined notch 


We now consider an inclined initial notch in the calculation domain in Fig. 4. The notch has an 
length of 0.8 m and an inclination angle 0 = 45° while the other simulation settings are the same 
as those in Subsection 4.1. Stress contrast o,9/ayo = [1, 2,4, 6, 8, 10] is applied with the vertical 
stress ayo = 0.5 MPa unchanged. Therefore, there are six simulations performed to show the 
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effect of stress contrast on the fracture propagation from the inclined notch. 

By using the proposed PFM, the phase field distribution for different o,9/ay0 at time t = 6 
s is shown in Fig. 14. As observed, the fracture propagates straight along the direction of the 
initial notch when oz0/o,) = 1 at time t = 6 s; however, if the stress ratio is larger than 1, 
the fracture deflects from the direction of the notch (0 = 45°). All the simulations consistently 
indicate that the larger the ratio of a0/oy0 is, the smaller angle the deflected fracture has from 
the direction of the maximum in-situ stress Smar (Czo) in this example. The fluid pressure-time 
curve for different 0,9/o,9 is shown in Fig. 15. Because the minimum in-situ stress oyo is kept 
constant, the fluid pressure at the center of the notch only increases slightly with the increasing 


stress contrast 0,0 / Oyo. 


4.3 Fracture from two perpendicularly crossed notches 


We also test the hydraulic fracture propagation from two perpendicularly crossed notches. The 
notches are located in the center of the calculation domain shown in Fig. 4 and both have a 
length of 0.8 m. The other simulation settings are the same as those in Subsection 4.1. We only 
consider four cases of in-situ stress as described in Table 3. This means that the direction of the 


maximum stress Smar and the minimum stress Smin varies in the four cases. 


Table 3: Four cases of boundary and in-situ stresses 


Case Horizontal stress Vertical stress 
Case 1 || czo = 1 MPa oo = 0.5 MPa 
Case 2 || ozo = 5 MPa 049 = 0.5 MPa 
Case 3 || czo = 0.5 MPa oyo = 1 MPa 
Case 4 || czo = 0.5 MPa oyo = 5 MPa 


By using the proposed PFM, the phase field distribution for different cases is shown in Fig. 
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a) Oz0/Fyo = =1 b) o20/oy0 = =2 c) 620/ey0 = =4 
d) ozo/Tyo = 6 e) Fx0/Tyo = 8 f) o20/ey0 = 10 


Fig. 14 Fracture propagation from an inclined initial notch at time t = 6 s under different 
o20/ Oyo 


Mid-point fluid pressure [MPa] 


t [s] 


Fig. 15 Effect of the ratio of o:9/oy9 on fluid pressure-time curve for an inclined initial notch 
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16. As observed, fractures only initiate and propagate from the notch perpendicular to the 
direction of the minimum stress Smin while the other notches do not even grow. Note that the 
direction of these fractures is consistent with the direction of Smin. Therefore, Fig. 16 validates 
the engineering observation that the fractures perpendicular to the direction of the minimum 
in-situ stress will initiate and propagate more easily. The fluid pressure-time curves for different 
cases are shown in Fig. 17. The curves for Cases 1 and 3 and those for Case 2 and 4 are almost 


the same because the maximum and minimum in-situ stresses are identical. 


) Case 1 ) Case 2 
) Case 3 ) Case 4 


Fig. 16 Fracture propagation from two perpendicularly crossed notches at time t = 5 s 
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Mid-point fluid pressure [MPa] 


t [s] 


Fig. 17 Fluid pressure-time curve for two perpendicularly crossed notch 
4.4 Fracture from two horizontal notches 


In this 2D example, hydraulic fracture from two horizontal notches is presented to further validate 
the proposed PFM. The geometry and boundary stress condition of this example are shown in 
Fig. 18a where the horizontal and vertical remote stresses are both 0.5 MPa for simplicity. The 
parameters and numerical settings are the same as Subsection 4.1. The phase field at = 2.5 
s is shown in Fig. 18b. A symmetrical fracture pattern is observed and the fracture deflection 
occurs during fluid injection. This observation is in good agreement with the “stress shadowing” 


phenomenon in engineering practice [53]. 


4.5 Linearly varying stress 


In this final 2D example, we test the effect of linearly varying stress field on hydraulic fracture 
propagation. The geometry of this example is the same as that in Fig. 4 while the boundary 
stress condition is shown in 19. Note that in this example, the gravity is applied in the vertical 


direction, and linearly varying horizontal stress acts in the horizontal direction, which can be 
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0.5 MPa 0.5 MPa 
| | | fos MPa 
(a) Geometry and boundary condition (b) Phase field at t = 2.5 s 


Fig. 18 Fracture propagation from two horizontal notches 


considered as a combination of self-weight stress field and tectonic stress field in underground 


geological environment. The parameters and numerical settings are the same as Subsection 4.1. 


2 MPa 


1 MPa 1 MPa 


| Gravity 


0.8 of notch 


5 MPa 5 MPa 


2.25 MPa 


Fig. 19 Linearly varying stress boundary condition 


The phase field evolution at different time is shown in Fig. 20. An asymmetric fracture 
pattern is observed and the fracture propagation is much easier towards the upper boundary 
than towards the bottom. This finding can be also verified in Fig. 21, which depicts the fracture 
length increment at different time. Figure 21 indicates that the fracture length from the upper 


tip of the pre-existing notch is much larger than that from the lower tip. The downwards 
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hydraulic fracture is hampered and even cannot propagate after t = 9 s. Figures 20 and 21 
strongly verifies that the hydraulic fracture tends to propagate towards the region with lowest 


fracture resistance. 


(a)t=4s (b) t=9s (c) t=14s 


Fig. 20 Fracture propagation under varying vertical stress field 
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Fig. 21 Fracture length increment under nearly varying stress field 


In summary, the 2D examples in Subsections 4.1 to 4.5 indicate the sensitiveness of the 
hydraulic fracture propagation to the stress boundary condition. Our proposed model, which 


involves the effect of initial stress field, is feasible and practicable in capturing the effect of 
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remote stresses on hydraulic fracture propagation and in producing the correct displacement 


field. 


5 3D example 


In this section, we test the performance of our method in modeling 3D hydraulic fracture prop- 
agation. Here, the last example is a 3D isotropic medium with a penny-shaped initial notch. 
The fluid source in the notch is set as 10 kg/(m?-s). The calculation domain is a cube with a 
dimension of 10 x 10 x 10 m? and has the same center as the initial notch, the height of which 
is in the z direction. The notch is parallel to the z-y plane and the radius and height are 0.8 m 
and 0.4 m, respectively. 

The parameters for calculation are identical to those in Table 2 except the length scale 
parameter is lọ = 0.25 m and the permeability ky is 5.21 x 107? m?. We employ 6-node prism 
elements to discretize the 3D domain while the maximum element size is set as 0.18 m for 
reducing computational cost. In addition, the time increment is set as 0.05 s in each simulation. 
In this 3D example, we only test the influence of the stress in the z direction 0,9. Therefore, ozo 
is set as 0.5, 1, and 5 MPa, respectively, while the stresses in the other two directions are fixed 
to 1 MPa. 

Hydraulic fracture propagation patterns in the 3D isotropic medium are shown in Fig. 22. 
It can be seen from Figs. 22a-f that for gz) = 0.5 MPa and 1 MPa, the fractures initiate and 
propagate only in the x-y plane while the area of the fractured domain increases as 0,9 decreases. 
Figures 22g-i indicate that when a, 9 is too large, the fracture propagation in the x-y plane is 


hindered and fractures only propagate in the z-z or y-z plane. 
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The fluid pressure-time curves under different czo are shown in Fig. 23. The fluid pressure 
is observed to drop after fracture initiation and the maximum fluid pressure increases with the 
increase in the in-situ stress 0,9. Comparing Figs. 22 and 23 indicate that the effect of in-situ 
stress on the 3D porous medium can be well captured by the proposed PFM in a fixed FE mesh 
without requiring any re-meshing or adaptive techniques. The 3D simulations fully verify the 
strong capability of our proposed PFM for predicting complex hydraulic fracture propagation in 


porous media subjected to stress boundary condition. 


32 


chinaXiv:202309.00096v1 


020 = 5 MPa 


= 
ee 


e ae > 
Ses IER Sho 
SNES NES 


(g)t=15s (h) t=21s (i) t=27s 


E ST 


<=> me 


orous medium with an interior notch under 


Mid-point fluid pressure [MPa] 


t [s] 


Fig. 23 Fluid pressure-time curve for 3D porous medium 


6 Conclusions 


A new phase field model for simulating quasi-static hydraulic fracture propagation in porous me- 
dia subjected to stress boundary conditions is proposed. A new energy functional is established 
to consider the effect of initial in-situ stress field. This energy functional is then used to achieve 
the governing equations for the displacement and phase fields through the variational approach. 
Biot poroelasticity theory is used to couple the fluid pressure field and the displacement field 
while the phase field is used for determining the fluid properties from the intact domain to the 
fully broken domain. 

The presented numerical examples in this work verify the capability of the proposed PFM in 
capturing complex hydraulic fracture growth patterns in 2D and 3D. The numerical examples 
also indicate that under stress boundary conditions the proposed approach can obtain correct 
displacement distribution and reflect the sensitiveness of hydraulic fracture propagation to stress 
boundary conditions. In future research, the proposed PFM can be more widely applied in HF 


practices where the stress boundaries are dominated and can be used to investigate the influence 
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of naturally-layered porous media and multi-zone HF on fracture propagation. 
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Appendix Finite element discretization 


We first derive the weak forms of all the governing equations as 


[ V+ aoo- apt) :e)a0= [fy suas (29) 


Qn 


f —2(1 — k)H(1 — 6)6¢dQ + f E (uvo -Voo+ 03o) dQ = 0 (30) 
Q Q lo 


f osian- | pv-vepan= f mas f (am — paxa) d (31) 
o ót Q aN Q Ot 


In an element with n nodes, the nodal values for the three fields (u, ¢, and p) are defined as 


ui, Qi, and p; (i = 1,2,--- ,n). The fields are then discretized as follows, 


u = X Niu;, g= > Ndi: p= De Nipi (32) 


where N; is the shape function at node 7. We then derive the gradients of the three fields as 
e=) Biu, Vo= > Bed, Vp= >_ Bep (33) 
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where B¥, B’, and B? are derivatives of the shape functions: 


T 
Na 0 O My Ü Ño Nin 
By=| 0 Ny O Na Nz 0 | > BP=BP=|N,, (34) 
0 0 Niz 0 N; y N; zr Niz 


For 2D, the components along the z direction is removed from above equations (34). Due 
to the arbitrariness in the test functions, the external force FY" and inner force FY" for the 


displacement field are described by 


Fie [vias + f Byapran — f BH glooo 
Qn Q S (35) 
ps f [B;]"od9Q 
Q 
The inner force term of the phase field is also obtained by 
; 1 
Forint = f —2(1 — k)(1 - ¢)HN; + G. (uBitvo + ZON: ) dQ (36) 
5 0 


f i ME i 
Finally, for the pressure field, the inner force FP”, viscous force F?’"’*, and external force 
2 ? a ? a ) 
t : 
FP" are given by 


K. 
Fp = f pat’ vpao 
Q Meff 


Bevis = f Nps a0 (37) 
Q 


o vo 
P f N; 2 — PAXR z ) dQ + NiMyds 
Q dQN 


Thus, contribution of node 7 to the residual of the discrete equations for the three field is 
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written as 


H = — Fein (38) 


j 


RP = peext = per = pees 


(3 7 7 7 


Because the staggered scheme is used to solve the displacement, phase field and fluid pressure 
sequentially. We also adopt the Newton-Raphson approach sequentially to achieve R¥ = 0, 


R? = 0, and R? = 0 for the three fields. In addition, the tangents on the element level are 


calculated by 


ORY / = 
Kw = —i_ = | [BY)™D,[B“]dQ 
J Ou; al ] | 4] 
Oren o 
KY = 3g, = [ {BACB +N; (21 Se +) nN} dQ (39) 
J 
on pK 
Ki = = fer TTB dQ 
Pj Q Leff 


where D, is the elasticity matrices derived from the elasticity tensor D = 00°/0e. 
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